/* ************************************************************************* **
**    OpenVFIFE - Open System for Vector Form Instrinsic                     **
**                Finite Element Method (VFIFE)                              **
**                GinkGo(Tan Biao)                                           **
**                                                                           **
**                                                                           **
** (C) Copyright 2021, The GinkGo(Tan Biao). All Rights Reserved.            **
**                                                                           **
** Commercial use of this program without express permission of              **
** GinkGo(Tan Biao), is strictly prohibited.  See                            **
** file 'COPYRIGHT'  in main directory for information on usage and          **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.                  **
**                                                                           **
** Developed by:                                                             **
**      Tan Biao (ginkgoltd@outlook.com)                                     **
**                                                                           **
** ************************************************************************* */

// $Date: 2020-05-11 $
// Written: Tan Biao
// Revised:
//
// Purpose: This file contains the class definition for materials, including
// base class BaseMaterial (abstract base class), UniaxialPlastic, NdPlastic,
// and several  derived materials - LinearElastic, NonlinearElastic,
// UniBilinearStatic, UniBilinear, UniPower.

// The interface:
//
// Caution: for nonlinear material, there will be numercial error which can be
// reduce by increase the increment of strain.


/*****************************************************************************/
/*                              Material Libary                              */
/* based on "Elasticity and Plasticity", chen & A.F. Salip                   */
/*****************************************************************************/
#ifndef MATERIAL_H_
#define MATERIAL_H_
#include <iostream>
#include <string>
#include <cmath>
#include "Eigen/Dense"


class BaseMaterial
{
    // when E, G and nu are all defined; the G = E / 2 / (1 + nu);
    protected:
        int id_;
        std::string type_;
        double E_;
        double G_;
        double nu_;
        double density_;
        double ce_;
        double damp_ratio_;
        double damp_coef_;
        double mu_;
        double yield_stress_;
        double limited_tensile_strain_;

    public:
        BaseMaterial(int id);
        virtual ~BaseMaterial();

        int id() const;
        std::string type() const;
        double E() const;
        double G() const;
        double nu() const;
        double density() const;
        double ce() const;
        double damp_ratio() const;
        double damp_coef() const;
        double mu() const;
        double yield_stress() const;
        double limited_tensile_strain() const;
        virtual bool isFractured(double strain) const;

        void setE(double val);
        void setG(double val);
        void setNu(double val);
        void setDensity(double val);
        void setDampRatio(double val);
        void setDampCoef(double val);
        void setMu(double val);
        void setLimitedTensileStrain(double tensile);
        virtual void info() const;

        double calcCe() const;
        /* parameters of Eq()
        * ds - delta strain;
        * stress - current stress,
        * strain - current strain,
        * ps - current plastic strain
        * alpha - back stress,
        * q - hardening funciton,
        */
        virtual double Eq(double ds=0, double* stress=nullptr,
            double* strain=nullptr, double* ps=nullptr, double* alpha=nullptr,
            double* q=nullptr) const = 0;
};


class LinearElastic: public BaseMaterial
{
    public:
        LinearElastic(int id);
        virtual ~LinearElastic();
        // will change stress and strain
        virtual double Eq(double ds=0, double* stress=nullptr,
            double* strain=nullptr, double* ps=nullptr, double* alpha=nullptr,
            double* q=nullptr) const;
};


class NonlinearElastic: public BaseMaterial
{
    public:
        NonlinearElastic(int id);
        virtual ~NonlinearElastic();
        virtual double constitutive(double a, ...) = 0;
        virtual double Eq(double ds=0, double* stress=nullptr,
            double* strain=nullptr, double* ps=nullptr, double* alpha=nullptr,
            double* q=nullptr) const = 0;
};


class UniaxialPlastic: public BaseMaterial
{
    // pure virtual class
    protected:
        double m_;  // hardening, 0 - kinematic, 1 - isotropic, 0~1 - mixed

        // determine if the material yield
        virtual bool isYield(double stress, const double* alpha,
                             const double* k) const;

    public:
        UniaxialPlastic(int id);
        virtual ~UniaxialPlastic();
        void setYieldStress(double val);
        virtual void info() const;
        virtual double Eq(double ds=0, double* stress=nullptr,
            double* strain=nullptr, double* ps=nullptr, double* alpha=nullptr,
            double* q=nullptr) const = 0;
};


class UniIdeal: public UniaxialPlastic
{
    protected:
        double Et_;  // Tangent Modulus, Et_ == 0
        virtual bool isYield(double stress, const double* alpha,
                             const double* k) const;

    public:
        UniIdeal(int id, double E, double yield_stress);
        virtual ~UniIdeal();
        // virtual void info() const;
        virtual double Eq(double ds=0, double* stress=nullptr,
            double* strain=nullptr, double* ps=nullptr, double* alpha=nullptr,
            double* q=nullptr) const;
};


class UniBilinear: public UniaxialPlastic
{
    protected:
        double Et_;  // Tagent Modulus

    public:
        UniBilinear(int id, double E, double Et, double yield_stress, double m);
        virtual ~UniBilinear();
        double Et() const;
        virtual void info() const;
        virtual double Eq(double ds=0, double* stress=nullptr,
            double* strain=nullptr, double* ps=nullptr, double* alpha=nullptr,
            double* q=nullptr) const;
};


class UniPower: public UniaxialPlastic
{
    // TODO: can't work! The Eq is not right(need to record last yield strain)
    /* stress = E * strain (not yield), stress = k * epsilon ** n (yield)*/
    /*
    * Equation: $\sigma = E \epsilon for \sigma <= {\sigma}_0$
    *           $\sigma = k {\epsilon}^n for \sigma > {\sigma}_0$
    *           ${\sigma}_0 = k({\sigma}_{0} / E) ^ n$
    */
    private:
        double k_;
        double n_;

    protected:

    public:
        UniPower(int id, double E, double n, double yield_stress, double m);
        virtual ~UniPower();
        virtual void info() const;
        virtual double Eq(double ds=0, double* stress=nullptr,
            double* strain=nullptr, double* ps=nullptr, double* alpha=nullptr,
            double* q=nullptr) const;
};


// class NdPlastic: public BaseMaterial
// {
//     // TODO: need to add more members
//     protected:
//         Eigen::Matrix3d Et_;
//         Eigen::Matrix3d Eq_;
//         Eigen::Matrix3d stress_;
//         Eigen::Matrix3d strain_;
//         Eigen::Matrix3d delta_strain_;

//         virtual void updateYieldFunc()=0;
//         virtual bool isYield()=0;
//         virtual bool isLoad()=0;

//     public:
//         NdPlastic(int id);
//         virtual ~NdPlastic();
//         virtual void info() const;
//         virtual Eigen::Matrix3d* Eq(Eigen::Matrix3d &delta_strain_) = 0;

//         virtual double Eq(double ds=0, double* stress=NULL, double* alpha=NULL,
//                   double* kappa=NULL, double* stress_history=NULL,
//                   double* strain_history=NULL) const = 0;
// };

#endif